function [ v , zeta , Lh , vU , converged ] =...
    SolveVL( Grids , Params , GE )

% Use bisection

ScaleFactor = Grids.ScaleFactorVLtop ;    
    
% Initialize bounds

vU0 = GE.vU ; 
vU1 = GE.vU ;

% Lower bound
GE.vU = vU0 ;
[ ~ , zeta , ~ ] = SolveODE( Grids , Params , GE ) ;
zetaBot0 = zeta(1) ;
enteredloop = 0 ;

while zetaBot0 > Params.zetaU
    vU0 = vU0 / ScaleFactor ;
    GE.vU = vU0 ;
    [ ~ , zeta , ~ ] = SolveODE( Grids , Params , GE ) ;
    zetaBot0 = zeta(1);
    enteredloop = 1 ;
end

% Upper bound
if enteredloop == 1
    GE.vU = ScaleFactor * vU0 ;
else
    GE.vU = vU1 ;
end
[ ~ , zeta , ~ ] = SolveODE( Grids , Params , GE ) ;
zetaBot1 = zeta(1) ;

while zetaBot1 < Params.zetaU
    vU1 = ScaleFactor * vU1 ;
    GE.vU = vU1 ;
    [ ~ , zeta , ~ ] = SolveODE( Grids , Params , GE ) ;
    zetaBot1 = zeta(1);
end

% Bisect

error = 1 ;
it = 1 ; 
while error > Grids.tolv && it < Grids.maxit
    
    vUmid = ( vU0 + vU1 ) / 2 ;
    GE.vU = vUmid ;
    [ ~ , zeta , ~ ] = SolveODE( Grids , Params , GE ) ;
    zetaBot = zeta(1) ;
    
    error1 = abs( Params.zetaU - zetaBot ) ;
    error2 = abs( ( Params.zetaU - zetaBot ) / Params.zetaU ) ;
    error = max(error1,error2) ;
    
    if zetaBot < Params.zetaU
        vU0 = vUmid ;
        zetaBot0 = zetaBot ;
    elseif zetaBot > Params.zetaU
        vU1 = vUmid ;
        zetaBot1 = zetaBot ;
    end
    
    it = it + 1 ;
        
end

frac = ( Params.zetaU - zetaBot0 ) / ( zetaBot1 - zetaBot0 ) ;
vU = ( 1 - frac ) * vU0 + frac * vU1 ;
GE.vU = vU ;
[ v , zeta , Lh ] = SolveODE( Grids , Params , GE ) ;

converged = it < Grids.maxit ;

end

